Incompatible pollen abortion and late-acting self-incompatibility in Schima superba

In angiosperms, self-incompatibility (SI) is a common and widespread mechanism for plant prevention of inbreeding, and late-acting self-incompatibility (LSI) may be ancestral in the group. In this work, we studied Schima superba, a species in Theaceae that is a commercially important timer and fire-resistant tree, and revealed its LSI mechanism. Hormones, enzymes, transcriptomes, and proteins were compared between self-pollination (SP) and outcross pollination (OP) in the styles and ovaries from 0 to 120 h after pollination. The self-pollen tubes grew to the bottom of the style and entered the ovary within 48 h but failed to penetrate the ovule. Meanwhile, the hormone and peroxidase levels dramatically changed. Transcriptome and proteome analyses explored the molecular mechanisms of LSI and candidate genes related to LSI in S. superba. Overall, 586.71 million reads were obtained, and 79,642 (39.08%) unigenes were annotated. KEGG and GO analysis showed that there were 4531 differentially expressed genes (DEGs) and 82 differentially expressed proteins (DEPs) at 48 h in self- (SP) versus outcross pollination (OP). Among these, 160 DEGs and 33 DEPs were involved in pollen–pistil interactions. “Pollen–pistil interaction,” “signal recognition,” and “component of membrane” were downregulated in SP, whereas “cell wall and membrane biosynthetic process,” and “oxidoreductase activity” were upregulated. The DEGs involved with S-RNases and SCF during SP suggested that the LSI occurred at 48 h in the ovary and that the LSI in S. superba was under gametophyte control. Calcium ion increase and release, mitochondrial function loss, and ROS disruption further aggravated PCD progress and cell death. The LSI of S. superba, which happened 48 h after pollination, was a key time point. The incompatibility PT ceased growth in the ovary because of S-RNase recognition and PCD in this organ. This study highlights the LSI molecular mechanism in S. superba and provides a reference to other species in Theaceae.

At 48 h after pollination, growth-related hormones, such as auxin, cytokinin, and ABA 31 , were significantly downregulated, and after this time, the ABA levels increased and exceeded the normal level ( Fig. 3A-C). At 48 h after pollination, the plant hormone signal transduction genes and proteins underwent corresponding changes;   Table S4).
In total, 82 DEPs were identified using KEGG analysis and functional enrichment, and 30 pathways were enriched, 10 of which were significantly enriched in metabolism and biosynthesis (Supplementary Table S7).

Association analysis between DEGs and DEPs.
Association analysis were performed to clarify the consistency of DEG and DEP expression patterns. DEPs did not have a significant correlation with DEGs at 48 h in the ovaries of SP versus OP (y = 0.047x + 0.1854, R 2 = 0.0436, P = 0.179, Fig. 7A). Only c113206_g1| SS8413 and c104041_g1| SS2022 were co-detected from the DEP-DEG association analysis (Fig. 7B). Forty-one proteins and genes annotated in 3 GO functional enrichments were aggregated into four clusters, and only cluster 3 was upregulated in cells (Supplementary Figure S6A). A total of 24 proteins and genes categorized into 7 KEGG pathways were also aggregated into four clusters (Supplementary Figure S6B). These pathways were mainly in oxidation-reduction and basal metabolism.

Discussion
SI is an important mechanism that protects flowering plants, allowing them to overcome inbreeding depression and providing a high level of heterozygosity 32 . Pollen tubes grow from the stigma to the ovaries, and penetrate the ovules. This forms a zygote and successfully promotes fruit and seed formation. A series of complex signaling controls then occurs. Incompatible pollen can be killed at any stage of this process. Based on our previous 30 and present studies, the morphological, anatomical, and gene expression patterns have been comprehensively explored to unravel the complexity of LSI in S. superba. The time point of 48 h, as implicated in this study, was important for PT elongation to the ovary and penetration of the ovule. The transcriptome and proteome data showed that DEGs and DEPs were more enriched in the "oxidation-reduction process, " "ribonuclease T2 activity, " "cell wall or membrane biosynthetic progress, " and "transmembrane transport" in ovaries of SP relative to ovaries www.nature.com/scientificreports/ of OP, while "recognition of pollen, " "plant hormone signal transduction, " and "Glycolysis" were downregulated (Supplementary Figure S2 and Supplementary Table S8).
Incompatible PT growth and its physiological changes. Researchers have already partly revealed the SI mechanism in some species in Theaceae 15,22,26,27,33 . Seth 27 and Chen 33 found that the PTs of SP and OP of C. sinensis entered the ovary at 48 h; however, the IPT recognition system was expressed in the style. Zhang 22 showed that PT growth in the SP of C. sinensis was hindered at 24 h in the style. The PT of SP and OP of C. oleif-  www.nature.com/scientificreports/ era entered the ovary at 60 h, but the PT of SP failed to penetrate the ovule. Therefore, 24 h for the style's SI or the ovary's LSI were two features in Theaceae. In our previous research, we found that PT growth stopped and that the IPT recognition system occurred in the ovary 48 h after pollination (Fig. 2) 30 . From the phylogenetic position of Theaceae and the fossil of Schima, it was shown that the Schima possessed n = 18 symmetric karyotypes, had 5, not 3 locules, had no center axis, and confirmed that the Schimeae (≡Gordonieae) belonged to an ancient taxon and that this branch should be an early node on the phylogenetic tree in Theaceae 24,25,29,34,35 . We found the ovary's LSI in S. superba was not effective; during self-pollination, 11.2% still developed fruit 30 , and C. sinensis had only 1.1% 27 . This may be for an IPT abortion in the ovary, and this distance may not have been enough to encounter more lethal factors than that from the style to the ovule. Reactive oxygen species (ROS) are the reactive products of oxygen that have the potential to damage living cells, and they play a key role in diverse development stages, such as self-incompatibility during pollination to induce PCD in incompatible pollen 36 . In plants, enzymatic and non-enzymatic antioxidants (proline, carotenoids, alpha-tocopherol, glutathione, ascorbic acid, flavonoids, and carotenoids) act as ROS detoxicants. In our results, GST-related proteins, such as the lambda GSTs (GSTL| SS2206) and dehydroascorbate reductases (DHARs| SS9106), which belong to the outlying minor GST classes, lactoylglutathione lyase (GLX1|SS3305, SS4306, SS4307), glutathione (GSH| SS7412| SS8421| c112627_g1| c86003_g1| c90752_g1| c97760_g1), and SOD were all downregulated at 48 h in the ovary of SP. The reduction of detoxification factor levels in SP induces ROS accumulation and damage that can be lethal 35-38 . PT guidance and IPT recognition. PT guidance is divided into two processes: pre-ovular and ovular guidance, both of which occur in the ovary 37 . Receptor-like kinases (RLKs) mediate signaling pathways to control ovular guidance 38,39 . In our results, 16 G-type lectin S-receptor-like serine/threonine-protein kinases, which included 5 RLK1 (c126066_g2| c100047_g1| c100047_g3| c103929_g1| c126606_g1) and one RKS1 (c106281_ g1), in addition to one PT reception gene at rapid alkalinization factor (RALF), c86998_g1, were significantly downregulated at 48 h in the SP ovary.
S-RNase, which is exclusively expressed in pistil, is associated with female determinants in GSI and implicated in genetically identical pollen and the rejection of self-pollen, and it would trigger mitochondrial collapse and IPT's PCD 40 . In our study, three S-RNase genes, which contain the catalytic domain of T2-type ribonuclease (PF00445), c110815_g1, c115031_g1(RNS1) and c94788_g1(RNS3), were significantly upregulated at 48 h in the ovary of SP. Comparing these three genes with other S-RNases in Rosaceae, Plantaginaceae, Rutaceae, Solanaceae, and Rubiaceae (Supplementary Figure S7), c110815_g1 was found to have a close relationship with Prunus T2/Stype RNase. When comparing the deduced amino acid sequences of these three S-RNases with the S-RNases reported for Prunus species, we identified introns I, II, and III, five conserved regions, and a hypervariable region located between C2 and C3 (Supplementary Figure S8). The three S-RNases possessed the HVa (RHV), HVb, and RC4 regions, which are considered prime candidates for the S specificity-determining region, mainly in Rosaceous S-RNases 41 . Interestingly, c110815_g1 had three introns: the first was next to the Prunus-specific intron and located in C1; the second was located in the same position as in Prunus S-RNase in RHV; and the third was in RC4. We presumed that this gene could have a special function in self-incompatibility.
The male determinant that interacts with S-RNase degradation is the S-locus region, which contains the novel F-box protein named SLF/SFBs 42 . The SLF/SFBs interacted with cognate S-RNase and prevented non-self S-RNase catalytic activity. When non-self S-RNase enters the PT, it forms an SCF complex and targets the S-RNase for ubiquitination and degradation; PT then continues growth. Self S-RNase binds to the recognition domain, not the active domain, resulting in dysfunction of SCF for polyubiquitination of self S-RNases. Then, self-S-RNases are www.nature.com/scientificreports/ released to trigger subsequent PT growth inhibition events 43 . In our results, 6 SFB putative genes were identified: F-box protein SKIP23-like (c88773_g1|c41648_g1) and F-box LRR-repeat protein (c120445_g3| c103038_g1| c51202_g1| c110213_g1). Phylogenetic analyses of S. superba SFB-like genes together with Rosaceae, Solanaceae and Plantaginaceae SFB and SKP genes, shown in Supplementary Figure S9, revealed that the c51202_g1 gene clusters with the Prunus SFB gene. Interestingly, c51202_g1 had 5 polymorphic SNP loci in the F-box domain, tightly linked to the RHV region (second intron) of c110815_g1 (Fig. 8). However, this gene was expressed in all tissues analyzed here (Supplementary Figure S10), whereas the S-pollen gene(s) were mainly expressed in anthers/pollen only, such as c41648_g1 and c110213_g1; thus, the c51202_g1 detected in S. superba, similar to Prunus species, may function as a general inhibitor (GI), as in Prunus, and c41648_g1 and c110213_g1 may function as SLF 40,44,45 . Further experiments should be conducted to clarify these proteins' interactions and verify this presumption. qRT-PCR showed that c110815_g1 was mainly expressed in the ovary and increased in SP ovaries from 24 to 72 h, and c51202_g1, c41648_g1, and c110213_g1 were significantly upregulated at 48 h in the ovary of SP (Fig. 9). We presumed that these candidate genes were active in the self S-RNase's cytotoxicity at 48 h in the ovary, and LSI in S. superba could be gametophyte controlled.

Programmed cell death in IPT.
Apoptosis or PCD is a highly conserved mechanism that removes unwanted cells in eukaryotes 46 . At 48 h in the ovary's LSI in S. superba, the IPTs ceased growth and died. As mentioned above, the S-RNase induces apoptosis or PCD by direct degradation of ribosome RNA or indirect phosphorylation of the protein, destabilization of cytoskeletons, and release of cytochrome c (cyt c) 47,48 . Indirect progress related to Ca 2+ signaling is transduced by protein phosphorylation, and this signal transduction mechanism relates to mitogen-activated protein kinase (MAPK) cascades. The Ca 2+ concentration in the cytosol can also trigger downstream sensors to form stimulus-specific information. In S. superba, one MAPK gene c120765_ g5 was significantly upregulated, and the "Ca 2+ channel activity", "calcium ion transmembrane transporter activity, " and "calcium ion transport" pathway-related genes and proteins were also upregulated. The important prophase characteristic of PCD inhibits mitochondrial outer membrane permeabilization (MOMP), and the mitochondria system loses its function 49 . Then, some soluble proteins, such as cyt c, diffuse from the intermembrane space (IMS) to cytosol 50,51 . GO analysis of DEGs showed that 30 GO class 14 DEGs related to mitochondria, such as inner mitochondrial membrane protein complex, mitochondrial membrane organization, and mitochondrial proton-transporting ATP synthase complex, were downregulated, and none of these were upregulated at 48 h in the ovary of SP (Supplementary Table S9). This demonstrated that, at 48 h in the ovary of SP, the mitochondrial membrane system was dysfunctional, and PCD had already started. In addition, during apoptosis, Bcl-2 family proteins regulate soluble proteins (e.g., cyt c), which are released from mitochondria 48,52 . The Bcl-2-associated athanogene-like protein gene c122719_g2, which contained the BAG domain, was 2.5-fold more highly expressed at 48 h in the ovary of SP (Supplementary Table S10). Kang reported that AtBAG6, which contains the BAG domain, was overexpressed in Arabidopsis, and promoted mitochondrial fusion, causing cyt c release into the cytosol and inducing PCD [52][53][54] .
Cyt c leakage from the mitochondria into the cytosol is an initial marker for PCD 47,48 . In our study, the concentrations of cytochrome oxidase COX (SS204) and ATP synthase (SS2015) were 9.7-and 5.8-fold higher, respectively, at 48 h in the ovary of SP. COX is the enzyme that catalyzes the oxidation of reduced cyt c by molecular oxygen. A higher COX concentration at 48 h in the ovary of SP indicates that the plants have already started the self-protection mechanism due to higher levels of cyt c. If the COX concentration exceeds the critical www.nature.com/scientificreports/ value, the activity is limited by either the cyt c or ascorbate concentration, and PCD occurs 59 . These results are consistent with the results of the disordered oxidation-reduction system shown previously.

Conclusions
In this study, the analysis of PT growth characteristics, internal physiology changes, transcriptome, and proteome of SP and OP ovaries revealed the LSI mechanism in S. superba and identified several candidate genes. The time point 48 h after pollination in the ovary's LSI in S. superba was different from 24 h after pollination in style's SI in C. sinensis. High expression levels of 2 S-RNase genes and 3 SLF-related genes suggest that LSI in S. superba is under gametophyte control. Self S-RNase induced PCD of IPT could occur by indirect transduction of Ca 2+ signaling and release of cyt c. Mitochondrial fusion and function loss caused the PCD; in addition, oxidation-reduction system disorder promoted cell death. These results revealed the LSI molecular mechanism in S. superba and provided a reference for other plants in the Theaceae family.

Plant materials and pollination treatment. This study was conducted by the Forestry Genetic and
Breeding Lab at the Research Institute of Subtropical Forestry, CAF, State Forestry Administration, China (RISF-CAF). The State Forestry Administration is responsible for national parks and other protected areas. No specific permission was required for these locations/activities, as they were based on a non-destructive collection of plant material. The species is not endangered or protected, and the locations are not privately owned or protected by law. Two cultivars, 'JO59' and 'YX1' , grown in Lanxi nursery, Zhejiang province, China, were used in this study. These cultivars were identified, bred, and preserved by our institution. The scions were collected from selected trees of the natural forest of S. superba from Jianou and Youxi in Fujian; they were then grafted on local stock and maintained in Lanxi until 2013. The collection of plant material complied with institutional, national, and international guidelines. Field studies were conducted in accordance with local legislation.
Then, the styles and ovaries of each sample with three replications were frozen in liquid nitrogen and stored at − 80 °C for enzyme activity, hormone, amino acid, protein, and RNA-seq analysis.
Enzyme assays. Three defense-related enzymes, peroxidase (POD), catalase (CAT), and superoxide dismutase (SOD), were analyzed using SP and OP ovaries from 24 to 120 h. POD activity was determined using the guaiacol-catalase activity method 55 . SOD activity was measured using inhibition in the photoreduction of www.nature.com/scientificreports/ nitro blue tetrazolium (NBT) 56 . CAT activity was measured spectrophotometrically by monitoring the decrease in absorbance at 240 nm 57 .
Hormone extraction and determination. Using HPLC-MS/MS, the levels of auxin (IAA), abscisic acid (ABA), zeatin (ZT) were determined by Zoonbio Biotechnology Co., Ltd. (Nanjing, China) 58 . The SP and OP ovaries (0.1 g fresh weight), with three biological replicates from 24 to 120 h, were selected to extract the hormone, and HPLC-MS/MS analysis was performed (Agilent 1290). HPLC analysis was carried out using a poroshell 120 SB-C18 column (2.1 mm × 150 mm; 2.7 μm). The parameters were set as follows: spray voltage, + 4500 V; automizing temperature, 400 °C; the pressure of the air curtain, nebulizer, and aux gas were 15, 65, and 70 psi, respectively. SAS v8.0 (SAS Institute, Cary, NC, USA) was used to analyze the data.
Transcriptome analysis. The total RNA from the S48 and O48 ovaries was extracted using an RNAprep pure Plant Kit (Dingguo, Beijing, China). RNA quality and quantity were verified using 1% agarose gels and Agilent 2100 (Agilent Technologies, CA, USA). The mRNA was enriched and purified by magnetic beads with oligo T (dT), cleaved and synthesized first-strand cDNA using random hexamers. The double-strand cDNA libraries were then purified using the AMPure XP system (Beckman Coulter, Beverly, USA). Transcriptome sequencing was performed using the Illumina HiSeq™ 2000 sequencing platform in Novogene (Tianjin, China). The adaptor and low-quality sequences (Q < 20 or less than 35 bp) were removed in raw reads. Then the clean reads were assembled using Trinity assembler 59 . The functions of the unigenes were annotated using Blastx searches with 1e-5 against the protein databases. Blast2GO software 60 was used for gene Ontology (GO) term analysis. The plant Transcription Factor Database (PlnTFDB) (version 3.0) was used to determine the transcription factors (TFs) 61 . FPKM was used to estimate the expression of each gene 62 , and the DESeq R package 1.10.1 was used to estimate the differential expression between the two treatments 63 . DEGs were determined with an adjusted p value < 0.05 determined by DESeq. The R-seq data were uploaded to the sequence read archive (accession no. SUB6596208).

Protein extraction and expression analysis.
With some modification, protein extraction was performed following Isaacson 64 . The SP and OP ovaries at 48 h were frozen and ground to power, and 1 g of the sample was used for protein extraction 65 . Then, 10 μg samples were run on 12% SDS-PAGE gel and visualized by CBB stain 66 . Gel images with 300 dots per inch were scanned using an image scanner (GE Healthcare, USA). A total of 450 μl solution with 1500 μg protein sample was used to finish isoelectric focusing with the following parameters: 50uA per strip, rehydration at 50 V for 8 h, 100 V for 1 h, 200 V for 1 h, 500 V for 1 h, 1000 V for 1 h, 1000-10,000 V (gradient) for 1 h, 10,000 V for 13 h, 500 V for 12 h, temperature, 20 °C. The Ettan-DALT-Six system was run for 45 min at 100 V and then at 300 V for 6-8 h. Using PDquest 8.0 software, all gel images were processed in three steps: spot detection, volumetric quantification, and matching. The differential protein spots were selected using two thresholds (p ≤ 0.05, fold change ≥ 2 or ≤ 0.5). An ABI 5800 MALDI-TOF/TOF Plus mass spectrometer (Applied Biosystems, Foster City, USA) was used for peptide MS and MS/MS detection, and Cal-Mix5 was used to calibrate the instrument (ABI5800 Calibration Mixture). GPS Explorer V3.6 software (Applied Biosystems, USA) with default parameters was used for data integration, and the proteins were identified using the MASCOT V2.3 search engine (Matrix Science Ltd., London, U.K.).
Transcriptome and proteome association analysis. The corresponding transcripts were identified using the gene ID of the proteome, and the corresponding relationship between the protein and the transcript was determined. A series of association analyses were then perfomed. Correlation analysis was conducted on the fold change (taken as log2) of genes (proteins) identified by transcriptome and proteome analyses in the two omics stydies. The drawing process was implemented using Python. The different multiples of the corresponding transcriptome genes were determined in the two omics studies. The difference multiples were taken as log2 and then plotted using the ComplexHeatmap of the R package. GO annotation of the transcriptome and proteome was used for GO function enrichment using the tool wego (http:// wego. genom ics. org. cn/). Clustering heat map analysis of GO/KEGG function enrichment was conducted using the fold change (taken as log2) of DEPs identified using the proteome and corresponding transcriptome genes. The drawing process was implemented using ComplexHeatmap. Ethics approval and consent to participate. All field studies were performed in accordance with the local legislation in China and complied with the convention on trade in endangered species.